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Abstract 

We examine numerical rounding errors of some deterministic solvers for systems 
of ordinary differential equations (ODEs). We show that the accumulation of 
rounding errors results in a solution that is inherently random and we obtain the 
theoretical distribution of the trajectory as a function of time, the step size and the 
numerical precision of the computer. We consider, in particular, systems which 
amplify the effect of the rounding errors so that over long time periods the solutions 
exhibit divergent behaviour. By performing multiple repetitions with different 
values of the time step size, we observe numerically the random distributions 
predicted theoretically. We mainly focus on the explicit Euler and RK4 methods 
but also briefly consider more complex algorithms such as the implicit solvers 
VODE and RADAU5. 

Abbreviated Title: Rounding errors in the numerical solution of ODEs 


1 Introduction 

Consider ordinary differential equations (ODEs) of the form 

x t = b(x t ). 


These can be solved numerically using iteration methods of the type 

x t+h = x t +f3(h,x t ) 


where f3(h, x)/h —> b{x) as h —> 0. 

The simplest example of this is the Euler method where /3(h, x) = hb(x). This 
method is generally not used in practice as it is relatively inaccurate and unstable 
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compared to other methods. However more practical methods, such as the fourth order 
Runge-Kutta formula (RK4), fall into this scheme. 

When solving an ordinary differential equation numerically, each time an iteration 
is performed, an error e is incurred due to rounding i.e. 

Xb h = X? + /3(h,X*)+e. (1) 

Rounding errors in numerical computations are an inevitable consequence of finite 
precision arithmetic. The first work thoroughly analysing the effects of rounding errors 
on numerical algorithms is the classical textbook jl^| . A recent comprehensive treatment 
of the behaviour of numerical algorithms in finite precision including an extensive list of 
references can be found in Q. Although rounding errors are not random in the sense 
that the exact error incurred in any given calculation is fully determined (see 0 or B), 
probabilistic models have been shown to adequately describe their behaviour. In fact, 
statistical analysis of rounding errors can be traced back to one of the first works on 
rounding error analysis ji]. 

Henrici jfjj; §] proposes a probabilistic model for the individual rounding errors, 
whereby they are independent and uniform, the exact distribution depending on the 
specific finite precision arithmetic being used. Using the central limit theorem he shows 
that the theoretical distribution of the error accumulated after a fixed number of steps in 
the numerical solution of an ODE is asymptotically normal with variance proportional 
to /r _1 . By varying the initial conditions he obtains numerical distributions for the 
accumulated errors, with good agreement. Hull and Swenson |TlJ test the validity of the 
above model by adding a randomly generated error with the same distribution at each 
stage of the calculation, and comparing the distribution of the accumulated errors with 
those obtained purely by rounding. They observe that although rounding is neither a 
random process nor are successive errors independent, probabilistic models appear to 
provide a good description of what actually happens. 

We shall concentrate on floating point arithmetic, as used by modern computers. 
However, our methods can be used equally well for any finite precision arithmetic. 
We use the model, discussed and tested by the authors above, whereby under generic 
conditions, the errors in m can be viewed as independent, zero mean, uniform random 
variables e* ~ U[— \Xj* i \2~ p , |X^ i |2 -p ], p being a constant determined by the precision of 
the computer. 

In the first half of the paper we analyse the cumulative effect of these rounding errors 
as the step size h —> 0. Where previous authors have considered the accumulated error 
at a particular point, we derive a theoretical model for the entire trajectory. Cases 
in M 2 where the ordinary differential equation has a saddle fixed point at the origin 
exhibit the most interesting behaviour, as the structure of the ODE system amplifies 
the effect of the rounding errors and causes the numerical solution to diverge from the 
actual solution. We show that in this case the solution X]* is inherently random and 
we obtain its theoretical distribution as an explicit function of time, the step size and 
the precision of the computer. We shall see that as the step size h —» 0, the numerical 
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solution exhibits three types of behaviour, depending on the time. More precisely, there 
exists a constant c, determined by the ODE system, such that for times much smaller 
than —cAogh the numerical solution converges to the actual solution; for times close 
to —clog/r the solution undergoes a transition, driven by a Gaussian random variable 
whose distribution we shall obtain; for times much larger than —clog/r the numerical 
solution diverges from the actual solution. 

In the second half of the paper, we perform numerical simulations which illustrate this 
behaviour. By performing multiple repetitions with different values of the time step size, 
we observe the random distributions predicted theoretically. Where previous authors 
have obtained their numerical distributions by varying the initial conditions, we do so 
by introducing small variations in the step size h. We show that during the transition 
period described in the previous paragraph, the numerical solution intersects straight 
lines through the origin and we compare the theoretical and numerical distributions for 
the points at which these intersections occur. Both the mean and the standard deviation 
of these distributions are of the form ah 7 , where 7 G (0,1/2] is a constant determined 
by the ODE system, and a can be found explicitly in terms of the precision of the 
computer, i.e. the number of bits used internally by the computer to represent floating 
point numbers. We mainly focus on the explicit Euler and RK4 methods, but show that 
the same behaviour is also observable for more complex algorithms such as the implicit 
solvers VODE and RADAU5. 


2 Theoretical background 

In the paper 0 , limiting results are established for sequences of Markov processes that 
approximate solutions of ordinary differential equations with saddle fixed points. We 
shall outline these results and then show that the rounding errors accumulated when 
performing numerical schemes for solving ordinary differential equations can be viewed 
as a special case of this. This enables us to quantify how the rounding errors combine, 
and show that the resulting numerical solutions exhibit random behaviour, the exact 
distribution of which is obtained. 

2.1 Behaviour of stochastic jump processes 

We are interested in ordinary differential equations of the form 


x t = b(x t ). ( 2 ) 

We focus on M 2 in the case where the origin is a saddle fixed point of the system i.e. 
b(x t ) = Bx t + r(x t ), where B is a matrix with eigenvalues A, —/x, with A,p > 0 and 
r{x) = 0(|a;| 2 ) is twice continuously differentiable. This case is of particular interest 
as the structure of the system amplifies the effect of the rounding errors and causes 
the numerical solution to diverge from the actual solution over large times. Similar 
behaviour can be observed in higher dimensions where the matrix B has at least one 
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positive and one negative eigenvalue, although the corresponding quantitative analysis 
is much harder and we do not go into it here. 

The phase portrait of ((21) in the neighbourhood of the origin is shown in Figure In 



Figure 1: The phase portrait of an ordinary differential equation having a saddle fixed 
point at the origin (taken from ( 3 ]). 


particular, there exists some Xq ^ 0 such that 4> t (x o) —> 0 as t —> oo, where </> is the flow 
associated with the ordinary differential equation ©• The set of such x 0 is the stable 
manifold. There also exists some x^ such that 0^ 1 (x oo ) —> 0 as t —> oo. The set of such 
Xoc is the unstable manifold. 

Fix an xo in the stable manifold and consider sequences of Markov processes 
starting from Xo, which converge to the solution of O over compact time intervals. 
The processes are indexed so that the variance of the fluctuations of X f N is inversely 
proportional to N. If we allow the value of t to grow with N as a constant times log N, 
X f N deviates from the stable solution to a limit which is inherently random, before 
converging to an unstable solution (see Figure I2D- More precisely, we observe three 
different types of behaviour depending on the time scale: 


A. 


On compact time intervals, Xjr converges to the stable solution of (J2J), the fluc¬ 
tuations around this limit being of order N~z. The exact distribution of the 
fluctuations is asymptotically N~ 2 ^ where 7 1 is the solution to a linear stochastic 
differential equation, described in 
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B. Let v\ and v -2 be the unit eigenvectors of B corresponding to —ji and A respectively. 
There exists some x 0 7 ^ 0, depending only on xo, and a Gaussian random variable 
Zoo, such that if t lies in the interval [ R , A log N — f?], then 

X = xq e ^(vi + ei) + N ^Z 00 e xt (v2 + 62 ) 
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Figure 2: Diagram showing how the Markov process Xf deviates from the stable solu¬ 
tion (j) t (x o) for large values oft (taken from 114] ). 

where e*(t, N) —> 0 uniformly in t in probability as R, N oo. In other words, 
can be approximated by the solution to the linear ordinary differential equation 

Vt = By t (3) 

starting from the random point XoVi + N~% Z OQ V 2 - 

C. Provided Z^ ^ 0, on time intervals of a fixed length around logIV, X t N con¬ 

verges to one of the two unstable solutions of ©, each with probability 1/2, 
depending on the sign of Z^. 

2.2 Accumulation of rounding errors 

We can apply the above results to describe quantitatively how rounding errors accu¬ 
mulate when solving ordinary differential equations of the form a numerically. In 
particular we consider using iteration methods of the type 


x t+ h = x t + /3{h,x t ) 

where f3(h, x)/h —> b(x ) as h —> 0 e.g. the Euler method where /3(h,x ) = hb(x). 

Each time an iteration is performed, an error e = e(h, t) is incurred due to rounding, 
so we obtain a process (X^) tem iteratively by 

X t h +h = X t h + (3(h,X t h ) + e. (4) 
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Modern computers store real numbers by expressing them in binary as x = m2 n for some 
1 ^ \m\ < 2 and n e Z, and allocate a fixed number of bits to store the mantissa m and 
a (different) fixed number of bits to store the exponent n B When adding a smaller 
number to x , the size of the rounding error incurred is between 0 and 2 n .2~ p = 2^ loS2 
where p is the number of bits allocated to store the mantissa. Although it is possible to 
carry out the calculations below using the exact value of 2 L log2 NJ - P ; the calculations are 
greatly simplified by approximating it by |a;|2 -p . This results in the ‘effective’ value of 
p differing from the actual value of p by some number between 0 and 1. Under generic 
conditions, the errors e can therefore be viewed as independent, mean zero, uniform 
random variables with approximate distribution e 

i rsJ U[-\X*A2->,\Xm~*] (see g E; 
0]). The assumption that the e* are independent is violated in certain pathological cases, 
for example where there is a lot of symmetry in the components. However, in general it 
is a reasonable assumption. 


Although the above iterations are carried out at discrete time intervals, it is conve¬ 
nient to embed the processes in continuous time by performing the iterations at times 
of a Poisson process with rate /r _1 . As /3(h, x) does not depend on t, this does not affect 
the shape of the resulting trajectories. In this way we obtain Markov processes Xp that 
approximate the stable solution of © for small values of h. If, in addition, we assume 


that 


h -i fWh .£> 



0 


as h —> 0 (note that both Euler and Runge-Kutta satisfy this condition), then under 
the correspondence N ~ h -1 , we satisfy the conditions needed to apply the results in 
B Our numerical solution therefore exhibits the following random behaviour: 


A. For times of order much smaller than — log h, approximates the stable solution 
of ©, the fluctuations around this limit being of order h%. 

B. There exists some Xq ^ 0, depending only on x 0 , and a Gaussian random variable 
Zoo, such that if t lies in the interval [—c log h, — ^ log h + c log h\ for some c > 0 , 
then X\ l is asymptotic to 

x 0 e~ pt vi +h^Zooe M V2, (5) 

the solution to the linear ordinary differential equation © starting from the ran¬ 
dom point x 0 Vi + h^Z oc V 2 - 

C. Provided ^ 0, on time intervals around — ^ylog/r whose length is of much 
smaller order than — log h, Xj 1 approximates one of the two unstable solutions of 
0. each with probability |, depending on the sign of Z^. 


The random behaviour resulting from the accumulation of rounding errors is most 
noticeable on time intervals of fixed lengths around — 2 (A+ At ) l°g as f° r these values 

of t the two terms XQe~ pt and 2 ’ 0 C e At in © are of the same order. Over this time 
interval, the numerical solution undergoes a transition from converging to the actual 
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solution to diverging from it. During this transition, for each value of 6 G (0, 7t/2), 
Xp crosses one of the straight lines passing through 0 in the direction cos^Ui ± sin 9v 2 . 
These intersections are important as they indicate the onset of divergent behaviour. 
The distribution of the point at which X\ l intersects one of the lines in the direction 
cos 9v i ± sin dv 2 is asymptotic to 

/ Z 2 (a+ m ) |Zoo|1xo| ^ | tan 9 |(cos 9v\ ± sin 9v 2 ). (6) 

In Section [231 we show how to evaluate the variance of doing so explicitly in the 
linear case and obtaining bounds in the non-linear case. In Section |3] we verify these 
results by numerically obtaining the predicted distribution for hitting a line through the 
origin. 

2.3 Explicit calculation of the variance 

Suppose that we are using a numerical scheme that satisfies the above conditions to 
obtain a solution to the ordinary differential equation © starting from Xq for some Xq 
in the stable manifold. In the non-linear case we require that Xo is sufficiently close 
to the origin such that t(xo) is small. In general, for simplicity, we shall assume that 
|x 0 | ^ 1. 

We define the ffow (j) associated with this system by 

4>t{x) = b{4> t {x)), 4> 0 {x) = 0 


and let x t = (pt{x o). 

Suppose V\ , v 2 G R 2 are the unit right-eigenvectors of B corresponding to — //, A 
respectively, and that v'^v^ G (R 2 )* are the corresponding left-eigenvectors (i.e. v[vj = 

*«)■ 

Define 

x 0 = lim e^v'^tixo) 


and 

D s = lim e“ A h4V0 t (x s ). 

t—>oo 

It is shown in [Q| that these limits exist and that |xq| ^ 21 xq | ^ 2, and \D S \ ^ 2. 


Finally, let 


a(x) 




be the covariance matrix of the multivariate uniform random variable e, defined in 
equation when X£ = x. Then ~ iV(0, cr^), where 


— 


e~ 2Xs D s a{x s )D*ds. 
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Note that o 2 ^ ^ ^2 2p . 

In the general non-linear case, evaluating cr^, explicitly is not possible as it involves 
solving 0- It is possible to obtain a better approximation than that above, although 
the important observation is that it is proportional to 2~ 2p . 

In the linear case, <j>t(x) = e Bt x and xo = |xoK ; i- Hence x t = |xo|e -M Vi, xo = |xo|, 
and D s = v' 2 , and so 

Note that the directions of V\ and v' 2 , relative to the standard basis, are critical. For 
example, if either v\ or v' 2 is parallel to one of the standard basis vectors, then = 0. 


3 Numerical experiments 

In this section we solve ODEs numerically using deterministic solvers and observe the 
predicted random distributions arising as a consequence of the accumulation of rounding 
errors. For simplicity, and in order to observe the desired effects as clearly as possible, 
we mainly focus on the most elementary of all numerical ODE solution methods: the 
standard explicit Euler algorithm with constant time step size. However, we observe 
similar behaviour for RK4 and also briefly mention results obtained with more com¬ 
plex solvers, such as VODE [l|. For a recent overview of ODE solvers see |2|; for an 
introductory text see 0]. 


3.1 The system 

For x : [0, oo) —> M 2 , consider the linear ODE 

x[t) = Bx(t) 


where 


B 




for fixed A , fi > 0. We introduce new coordinates 


x(t) = R(ip)x(t) 

by rotating about the origin by a fixed angle tp 6 [0, vr/2), i.e. 


R(<p) 


f cos tp — sin <p 
^sint^ cos p 


We arrive at the transformed system 


x{t) = B(p)x(t) 


( 7 ) 



with 


B(<p) = R{<p)BR(ip) T , 

which will be the system under consideration in the following 
initial value 

,(0)=*(,) (i) = (zz ). 

The phase space evolution is sketched in Figure 01 



Throughout, we use as 
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Figure 3: Phase space for the saddlepoint ODE system © with sample trajectories and 
lines where hitting distributions are recorded (dashed lines). 


3.2 Theoretical hitting distribution 

As discussed in Section E3 the numerical solution to the above ODE system undergoes 
a transition from converging to the actual solution to diverging from it. During this 
transition, the numerical trajectory crosses one of the straight lines passing through 0 
at an angle </> ± 9 for each value of 6 e (0,7r/2). These intersections are important as 
they indicate the onset of divergent behaviour. The hitting distributions also provide a 
means of measuring the random variable Z^, which drives the random variations in our 
solutions, and hence of verifying the theoretical results. 

Equation © gives the asymptotic distribution of the magnitude of the point at which 
the numerical solution hits the line through the origin at angle ip ± ^ as Z \ x +i‘ where 
Z is a Gaussian random variable with mean 0 and variance 

= h °L = 3 ^ 1 +n ) h2 ~ 2p ( COS( ^ sin ^) 2 ( 9 ) 
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i.e. Z ~ A/"(0,<t 2 ). We obtain an explicit formula for the asymptotic distribution by 
starting from the A/"(0, a 2 ) distribution 


p(x)dx = 


\J~2lTi 


exp l — 


7T(J 


2(X 2 ' 


x 2 Ida; 


and performing a transformation of the variable given by y — |a;| A +' j . The result is: 


( \ A 2(A + p) a 

p{y)ay = —7= —2/ M ex P 

V 2nap 


1 2(A+^) \ 

y 11 Jdy. 


2a 2 


In the case A = (i = 1 which we shall consider below, setting a = we obtain the 

family of distributions 


f(x)dx = ax ex p ^ — — a 2 x 4 ^jdx (10) 

which we shall fit to our numerical data to confirm our theoretical value of a. 


3.2.1 Choice of parameters 

Since there is no random element contained in a deterministic solver, for each repetition 
at least one parameter has to vary. At first, we tried varying the initial value in the 
direction of the eigenvectors, but this did not yield any interesting results. The chosen 
distribution of initial values was reproduced exactly in the hitting distribution and no 
randomness could be observed. Also, from an aesthetic point of view, it is preferable 
to vary only ‘internal’, i.e. numerical parameters of an algorithm, such as the time step 
size or error tolerances, instead of varying ‘physical’ parameters of the system like the 
initial value. 

The only internal parameter of Euler’s algorithm is the time step size h, which we 
varied as follows. Given a user-supplied value of h, the step size h, for the i th repetition 
(i G {!,..., L }) is defined by 


hi = h + A h(i — 1 — k), 

where the number of repetitions L = 2k + 1 and 0 < Ah -C h are also user-supplied. 
For all simulations, we set k = 10 4 . 

In each run, the trajectory at some point intersects the lines X\ = ±X2 (the dashed 
lines in Figure 01) • In order to produce histograms, we partition the interval [0,1] into 
a given fixed number of subintervals of equal length and count how many times y falls 
into each subinterval, where y denotes the distance of the point of intersection from the 
origin. 

Since the limit distribution is given by |Z| A +^, if the values of A and y differed 
significantly then the distribution would be hard to observe in a numerical experiment. 
Furthermore, if A /i then the trajectories are very quickly pushed away from the X\- 
axis so that fluctuations (i.e. rounding errors) become largely irrelevant. Conversely, if 
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//> A then the amplification of deviations from the aq-axis is too weak to be observable. 
These arguments suggest choosing A and /i within the same order of magnitude, and we 
therefore chose A = /i = 1 for all simulations. 

Another subtlety concerns the choice of the rotation angle ip . For certain values, 
trivial trajectories or symmetry effects can occur which conceal the desired accumulation 
of rounding errors. For instance, for p = 0 the second component X 2 of the solution is 
always zero, and therefore the trajectory stays on the line x 2 = 0 (or equivalently X 2 = 0) 
with no fluctuations. Note that this is in agreement with a 2 = 0 in equation (JHJ). 
For (p = 7t/ 4, any rounding error that appears in one component also appears in the 
other one, which implies that, again, the trajectory always stays on the line X 2 = 0 
(or equivalently X\ = x- 2 ). This case is pathological as it consistently violates our 
assumption that the rounding errors are independent for each component. For these 
reasons, we chose p = 7t/ 5 throughout. 

Reasonable choices of h and Ah are limited by several factors. If h is too large (in 
the considered case h > 1CF 1 for both single and double precision) then the observed 
hitting distributions differ substantially from the theoretical one because not enough 
rounding errors can accumulate and hence the deviations are not sufficiently random. 
The onset of such effects can be seen for large values of h in Figure EH Lower bounds 
on h are imposed on the one hand by computational cost and on the other hand by 
the numerical accuracy of the computer on which the calculations are performed. In 
practice, however, computational expense becomes prohibitive for values of h much 
larger than the smallest values permitted by numerical accuracy. Our particular choice 
of step size distribution requires that kAh should be (much) smaller than h. The lower 
limit for Ah is determined solely by the numerical precision, i.e. Ahjh must not be 
smaller than the numerical precision because otherwise there is no variation at all and 
hence no randomness. 

It is beyond the scope of this paper to investigate in detail the dependence of our 
observations on the distribution of step sizes. However, preliminary experiments with 
varying Ah and even with non-uniform step size distributions suggest that this depen¬ 
dence is very weak for a wide range of conditions. The mean value of the step size 
distribution, whose effect on the shape of the hitting distribution (i.e. the parameter a 
in equation (HOjl ) has been demonstrated in Figure O seems to be most significant. On 
the other hand, the variance seems irrelevant apart from being large enough to produce 
randomness, and does not affect the parameter a. This is also supported by the fact 
that a is (asymptotically) independent of the number of repetitions L. 

Figure 0] shows that the shape of the distribution exhibits no discernible systematic 
dependence on Ah over at least nine orders of magnitude. This reinforces that the 
parameter a is unaffected by the variance of the step size distribution for this particular 
case. The deviations seen for values of Ah smaller than about 10~ 19 are due to the fact 
that Ahjh approaches the limits of numerical precision. 
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Figure 4: Step size variation for Euler’s algorithm (double precision, step size h = 10 4 , 
L = 20001 repetitions each). 

3.2.2 Results and observations for explicit methods 

Using the given values, distributions as shown in Figure E] are obtained. We fitted the 
theoretical distribution m to the ones produced numerically with very good agreement. 

From the numerical experiments we obtain for each h a distribution of the form 
m, where the parameter a is given as a result of the fitting procedure. In Figure 
El the parameter a is plotted as a function of the time step size h, both for single 
(Figure |6( a) | ) and double (Figure |6(b)| ) precision (4 and 8 bytes internal representation 
of floating point numbers respectively). Error bars due to the fit are only about 1% 
and hence invisibly small. In both cases, the dependence between a and h seems to 
be well described by a oc \[Ji. Intuitively, one might explain the qualitative behaviour 
as follows: The smaller the step size, the more numerical errors accumulate, and hence 
the broader the distribution. The intuition “the smaller the step size, the higher the 
accuracy, and hence the narrower the distribution”, although at first sight equally valid, 
is false. In the single precision case, because of the lower accuracy compared to double 
precision, the distributions are much broader for a given step size. 

Equation m predicts the value of ah a to be 

= 4V ^ -x 2 P = 8.220 x 2 P . 

\pn cos | sin | 

For Euler’s method, the above data give ah~i = 9.411 x 10' for single precision and 
ah~ 2 = 4.956 x 10 16 for double precision. For 4 th order Runge-Kutta, the values are 
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Figure 5: Observed hitting distributions with theoretical fits for Euler’s algorithm (Ah = 
10 -10 , L = 20001 repetitions each). 

ah ~^ = 9.27 x 10' (with a relatively large error of ±0.12 x 10') for single precision and 
ah~ 2 = 4.746 x 10 16 for double precision. Using the approximation discussed in Section 
0 the actual value of p is between 23 and 24, when working in single precision, and 
between 52 and 53 when working in double precision, the particular value depending on 
the exact number being computed. Our theoretical results therefore predict ah~ 2 lies 
between 6.895 x 10' and 1.379 x 10 8 for single precision and between 3.702 x 10 16 and 
7.404 x 10 16 for double precision. 

There are three possible sources of error in our calculations. The first is the error 
in fitting the numerical data to the theoretical model, the second is that our theoretical 
models are based on asymptotic results as h —> 0, whereas we are applying them to 
values of h which are necessarily larger than the precision of the computer. The third 
source of error arises from the assumption that at each stage the rounding error can be 
viewed as an independent uniform random variable, depending on a fixed value of p. 
The above results show the above errors are all small and our theoretical model provides 
a very good fit. 

3.2.3 Implicit solvers 

Possible internal parameters to be varied in solver packages more sophisticated than 
Euler’s method are typically the error tolerances RT0L (relative) and AT0L (absolute) 
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(a) Single precision (Ah = 10 8 ). (b) Double precision (Ah = 10 10 ). 

Figure 6: Parameter a in equation (HUD as function of the time step size h for simple 
explicit methods (Euler and 4 th order Runge-Kutta). 


and the global time step h g (the time interval after which the user requests solution 
output from the solver), ffowever, naturally the user has no immediate control over 
the size of the actual steps taken, which is determined algorithmically as a function of 
the error tolerance parameters RTOL and ATOL, frequently by trial-and-error methods 
using heuristics, but only rarely by an explicit formula. Nonetheless, as shown in Figure 
[Tj distributions very similar to the ones seen for Euler’s algorithm (Figure 0) can be 
generated. Experiments do not readily suggest a simple relationship between the shape 
of the distribution (parameter a in equation fTTfl) ') and any of the parameters ATOL, RTOL, 
and h g . We suspect the lack of direct control over the time step size to be the main 
reason for this behaviour. We found that in order to produce Figure 0 one has to use 
RTOL = 0, which we also attribute to the step adaptation. 

For the solver RADAU5 |5] the results are qualitatively similar, which supports the 
assertion that the observed phenomena are not specific to a particular algorithm, but 
rather general effects. 


4 Conclusion 

We analysed the cumulative effect of rounding errors incurred by deterministic ODE 
solvers as the step size h —> 0. We considered in particular the interesting case where 
the ordinary differential equation has a saddle fixed point and showed that the numerical 
solution is inherently random and also obtained its theoretical distribution in terms of 
the time, step size and numerical precision. We showed that as the step size h —> 0, the 


14 



































Figure 7: Hitting distributions for VODE. 


numerical solution exhibits three types of behaviour, depending on the time: initially it 
converges to the actual solution, it then undergoes a transition stage, finally it diverges 
from the actual solution. 

By performing multiple repetitions with different values of the time step size, we 
observed the random distributions predicted theoretically. We demonstrated that during 
the transition period described above the numerical solution intersects all the straight 
lines through the origin. The theoretical and numerical distributions for the points 
at which these intersections occur showed very good agreement. Both the mean and 
the standard deviation of these distributions were found to be of the form a/r 7 , where 
76 (0,1/2] is a constant determined by the ODE system, and a was found explicitly in 
terms of the precision of the computer. We mainly focused on the explicit Euler and RK4 
methods, however, we also briefly considered the implicit solvers VODE and RADAU5 to 
demonstrate that the observed effects are not specific to a particular numerical method. 
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